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Abstract 


Estimates of the total particulate mass of the plumes of Enceladus are impor¬ 
tant to constrain theories of particle formation and transport at the surface 


and interior of the satellite. We revisit the calculations of Ingersoll and 


Ewald (2011), who estimated the particulate mass of the Enceladus plumes 


from strongly forward scattered light in Cassini ISS images. We model the 
plume as a combination of spherical particles and irregular aggregates result¬ 
ing from the coagulation of spherical monomers, the latter of which allows for 
plumes of lower particulate mass. Though a continuum of solutions are per¬ 
mitted by the model, the best fits to the ISS data consist either of low mass 
plumes composed entirely of small aggregates or high mass plumes composed 
of mostly spheres. The high particulate mass plumes have total particulate 


masses of (166 ± 42) x 10 3 kg, consistent with the results of Ingersoll and 


Ewald (2011). The low particulate mass plumes have masses of (25 ± 4) x 
10 3 kg, leading to a solid to vapor mass ratio of 0.07 ± 0.01 for the plume. If 
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indeed the plumes are made of such aggregates, then a vapor-based origin for 
the plume particles cannot be ruled out. Finally, we show that the residence 
time of the monomers inside the plume vents is sufficiently long for Brownian 
coagulation to form the aggregates before they are ejected to space. 
Keywords: Enceladus, ices, photometry, radiative transfer, satellites; 
surfaces 


1. Introduction 


Enceladus’ plumes provide an indirect way to study the subsurface. In 
particular, the ratio of ice particle to vapor mass can serve as an important 
constraint of ice particle formation and transport (Ingersoll and Ewald 20ll 


henceforth IE11). IE11 examined Cassini Imaging Science System (ISS) im¬ 
ages of the plumes at small scattering angles from 2.2 — 5.3° , where forward 
scattering is dominated by ice grains as opposed to water vapor. They then 
fit the resulting phase curves to various mass and shape distributions, assum¬ 
ing that the ice grains are solid spherical or ellipsoidal particles. However, 
the large particle to vapor mass ratio calculated by IE11 could not be easily 
explained by any existing theoretical models. Furthermore, recent results 
from the Cassini Cosmic Dust Analyzer (CDA) (Kempf, S., Cassini Project 
Science Meeting, Jan 22, 2015, and private communication) indicate that the 
particulate mass may be a factor of ten lower than the estimates of IE11. 
Since the defining property of the plume particles is strong forward scatter¬ 
ing, an alternative model for the plume is one made up of aggregates. 

An aggregate is a particle of irregular shape composed of smaller subunits 
(or “monomers”) stuck to each other. The monomers are usually composed 
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of simpler shapes, such a spheres, plates, or columns, depending on the for¬ 
mation mechanisms. Aggregates form under many different settings. For 
example, the hazes in the atmospheres of Jupiter and Titan are thought to 


be composed of hydrocarbon aggregates (West and Smith 1991 Tomasko 


et al. 2008 Zhang et al. 2013); Saturn’s F-ring is likely populated by ice 


aggregate particles (Vahidinia et al. 2011); and cirrus cloud particles on 
Earth take on a variety of non-spherical shapes, ranging from fernlike and 


fractal geometries to aggregates of irregular shapes (Yang and Liou 1998). 
In particular, ice clouds that form as a result of strong vertical motions are 
dominated by aggregate particles (Baum et al. 2011| Heymsheld et al., 2002). 
Since aggregates form under such a wide variety of conditions, it is plausible 
that they could form in the plumes of Enceladus. 

Thanks to the diversity of instruments on Cassini, the plumes are well 
studied. There are good estimates of water vapor mass from the Ultraviolet 
Imaging Spectrograph (UVIS) instrument (Tian et al., 2007; Hansen et al.| 


2011, for example), as well as that of other minor gaseous constituents, such 


as carbon dioxide, methane, ammonia, and argon from the Ion and Neutral 


Mass Spectrometer (INMS) (Waite Jr et al. 2009). Data from Cassini CDA 
suggests that the ice particles in the plumes can be broadly classified into 
two types: slow, large, salt-rich grains that tend to fall back onto the surface, 


and fast, salt-poor grains that escape into the E-ring (Postberg et al. 2011). 
Ice particle velocity distributions were measured using the Visual and In¬ 


frared Mapping Spectrometer (VIMS) (Hedrnan et ah, 2009J) and ISS (IE11). 
However, the shapes of the particles remain mostly unconstrained by obser¬ 
vations and are usually assumed to be spherical or oblate/prolate, such as in 
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Porco et al. (2006) and IE11. In this study, we derive estimates for the total 


particulate mass of the Enceladus plumes by extending the range of possible 
particle shapes to aggregates. 


2. Aggregate Model 

Aggregate particles are defined by two parameters: their fractal dimension 
D and their monomer radius r m . These two quantities are related to the 
number of monomers that make up the aggregate N m and the radius of the 
aggregate particle r by 

Nm=(^-) D . ( 1 ) 


A typical value of D for aggregates in the Solar System is around 2. For 
example, it has been shown that D — 2 is a good approximation for the 
aggregate particles in the Titan hydrocarbon haze, where the actual dimen¬ 


sion may vary between 1.75 and 2.5 (Cabane et al. 1993). D = 2 is also 


appropriate for snowflakes and cirrus cloud ice crystals on Earth, which have 


variations in D between 1.9 and 2.3 (Westbrook et al. 2006; Schmitt and 


Heymsheld, 2010). Zhang et al. (2013) further showed that D — 2 aggre 


gates can be used to fit Cassini ISS observations of the Jupiter stratospheric 
aerosols. Such particles have masses that scale linearly with surface area, like 
a sheet, though the particle itself is a three dimensional object. As a result, 
these particles tend to have small masses associated with large scattering 
cross sections. 

By comparison, r m values vary considerably across different types of ag¬ 
gregates. For example, r m = 10 nm for the stratospheric aerosols of Jupiter 


'Zhang et al. 

2013 

) and 40 nm for Titan’s haze aggregates ( 

Tomasko et al. 
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2009), but these values are easily dwarfed by that of Saturn’s F ring parti¬ 


cles, which can reach a few microns (Vahidinia et ah, 2011), while ice crystal 


monomers on Earth can be hundreds of microns across (Kajikawa and ffeyrns- 


held, 1989). We will therefore leave r m as a free parameter in our model that 


will be varied to best fit the data. As a simplification, we assume that all 
aggregates in the plume have monomers of the same r m . 

Equation [T| leads to a minimum size r min for an aggregate of 


r — 9 D T 
min ^ 1 m 


( 2 ) 


where we have chosen N m = 2 as the minimum number of monomers an 
aggregate can have. Particles with r < r min are assumed to be spherical with 
radius r. To simplify the problem and reduce the number of free parameters, 
we further assume that both spherical and aggregate particles “share” the 
same particle size distribution 


— = N 0 (r/r 0 ) f 3 /[l + (r/r 0 ) 2/ ], (3) 

which is the number of particles in the natural log of radius interval dln(r), 
with No as a parameter that scales with the particle number density, / as 
a positive width factor, and r 0 as the median radius of the particle mass 
distribution given by 


where 


dM(r) _ 2 Mo /(r/ro)^ 1 
dr 7ir 0 1 + [rjr o) 2 ^ ’ 


M (r) 


2 M 0 

7T 


arctan 



( 4 ) 

( 5 ) 
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is the total mass of particles with radius between 0 and r, and M 0 is the 
total mass of particles, ro splits the mass distribution into equal halves, with 
the width of the distribution governed by /; small / values indicate wide 
distributions while large / values indicate narrow distributions. M 0 , r 0 , and 
/ are free to vary in the model during optimization of the fit to the data. Eqs. 
[3]- [5] are the same distribution functions given in IE11, which were chosen due 
to their relative simplicity and their ability to capture both sharply peaked 
and asymptotic functional forms with only two free parameters. Figure [l] 
shows a schematic of how the aggregate and spherical particles “share” the 
dM{r)/dr distribution in our model. Particles with r > r min are assumed 
to be aggregates with radius r given in Eq. [l] and r m defined by Eq. [2j 
they follow the size distribution of Eq. [3] with some given ro and / values. 
Particles with r < r m i n are assumed to be spheres of radius r that follow the 
same size distribution as the aggregates, with the same r 0 and / values. With 
these definitions, dM{r)/dr has a discontinuity at r = r min . This is caused 
by the different ways the mass of single particles (aggregate or spheres) scales 
with r while keeping N 0 fixed for both the aggregate and spherical sections 
of the size distribution. The discontinuity is not shown in Figure [lj and we 
avoid it in our calculations, as discussed below. 

3. Observations and Model Setup 

Following the procedure of IE11, we use the following relationship to 
estimate total plume particulate mass (Eq. 4 of that paper): 



/ 


A p Q sca P(9)(dN/dln(r))dln(r) 
V p (dN / dln(r))dln(r) 


( 6 ) 
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where M 0 is now the total particulate mass of the plume; A p Q sca is the 
scattering cross section of the particle, which, for spherical particles, can 
be split into the geometric cross section of the particle A p = 7rr 2 and the 
scattering efficiency Q sca ; P(&) is the scattering phase function; 6 is the 
scattering angle, which is given in Table 1 of IE11; V p is the solid volume of 
each particle, given by (4/3)7rr 3 for spheres and (4/3)7rr^A^ m for aggregates; 
Pice = 0.917 g cm -3 is the density of ice; and r is the particle radius, as given 
for aggregates and spheres in Section [2j R{6 ) is defined as 


R{8) = j jiA 


( 7 ) 


where / is the measured radiance, nF is the solar irradiance, and the integral 
is taken over the area A of the image projected onto a plane at the distance 
of Enceladus to Cassini. The values and associated uncertainties of R(9) are 
also listed in Table 1 of IE11. Only the wide-angle camera (WAC) images are 
used. In order to remove the bias of the data towards scattering angles where 
more observations were taken, we repeat the procedure of IE11 and grouped 
the 18 brightness measurements into 10 data points by scattering angle, and 
averaged the brightness values and uncertainties within each group. 

In order to sample multiple scattering angles the observations were taken 


at multiple orbital phases. However, Hedman et al. (2013) showed that the 
brightness of the plumes varied as a function of orbital phase, even when 
scattering was taken into account, indicating variability in the intensity of 
the plumes. They showed that the plume brightness increased rapidly from 
an orbital phase of 90° to 180°, while below 90° the changes in brightness were 
much more gradual (see their Figure 4). Our observations, as given in Table 
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1 of IE11, were taken at orbital phases of 7.68° to 55.93°, well within the 
“gradual” part of the orbit. Therefore, while the intensity of the plume still 
changed during our observational period, we will assume that the changes 
are minor compared to brightness changes due to particle scattering. 

We calculate the Mie scattering cross section and phase function of the 


spherical particles using Grainger et al. (2004), while those of the aggregates 


are calculated using the model of Tomasko et al. (2008, henceforth T08), a 
detailed description of which can be found in the appendix of that paper. 
The T08 model is a fast parametrization valid only for D = 2 due to the 
computationally intensive calculations required for a model of aggregates of 
arbitrary dimension (Mishchenko et al. ,[l996); as discussed in Section [2j D 
= 2 is valid for many types of aggregates in the Solar System, including ice 
particles in Earth clouds. A visual description of such particles is found in 


Figure 1 of West and Smith (1991). The aggregate scattering model considers 
each aggregate to be a collection of spherical monomers that scatter and 
absorb as Mie spheres; the total intensity of the scattered light from the 
aggregate is then the superposed intensity of scattered light from all the 
monomers in the aggregate. For more complex interactions, such as repeated 
scattering and absorption by multiple monomers, correction factors are added 
that are functions of N m , the refractive index, and the size parameter of the 
monomers 




27Tr m 

A 


( 8 ) 


where A is the wavelength of the radiation being scattered. The correction 
factors are empirical and listed in Table A2 of Appendix A of T08. Only 














cases with x m between 10 -4 and 1.5 have currently been validated; and ex¬ 
trapolations to different x m values become increasingly uncertain the farther 
they are from the tested range of values (see again Appendix A of T08). The 
T08 model has been used to fit observations of Titan’s haze, as well as that 


of the stratospheric aerosols of Jupiter (Zhang et al. 2013) and aggregate 


particles created in laboratory experiments (Bar-Nun et ah, 2008). We as¬ 
sume a real index of refraction of 1.31 and an imaginary index of refraction 
of 0.0005 for ice in our scattering calculations for the three Cassini ISS wave¬ 
length channels considered (VIO, 420 nm; CLR, 635 nm; IR3, 918 nm), as 


with IE11 (Meng et al. 2010). 

In order to avoid the discontinuity in dM(r)/dr, we treat particles of 
different sizes separately by discretizing dN/dlnr into 25 bins in our calcula¬ 
tions, with a minimum particle radius of 80 nm and a particle volume ratio 
of 2 between successive bins. We then approximate the integrals in Eq. [6] as 
summations over the size bins, so that different forms of A pi Q sca , P(0), and 
V p can be used depending on whether the bin corresponds to spherical parti¬ 
cles or aggregate particles. This enables each integral to sum over both types 
of particles across the entire size range, provided N 0 is a constant that drops 
out of Eq. [6j and ensures that M 0 is a constant that applies to the entire 
plume, as Mq = pi ce f V p (dN / dln{r))dln{r) . The minimum particle size is set 
to prevent an infinite integral of dN/dlnr for / < 3 in our numerical scheme, 
though the total mass and R{9) are always finite for any positive value of /. 
This is because V p scales as r 3 for spheres and r 2 for aggregates, and because 
Qsca scales as r 4 when r —> 0. The value of the minimum particle radius was 
picked such that it is much lower than the likely range of favorable r m values 
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estimated by fitting the data by hand. We have also confirmed, through tests 
of our model, that particles smaller than 80 nm do not contribute signhcantly 
to the forward scattering for all reasonable model parameter values. 


4. Fitting to Data 


We use the Markov Chain Monte Carlo (MCMC) method (Hastings, 1970) 
to optimize the fits to the data and obtain statistically significant confidence 
intervals for the optimized parameters, which consist of M 0 , r m , r 0 , and /, 
with D fixed to 2. This is the same number of parameters as IE11, though 
they had a particle aspect ratio rather than r m . 

The MCMC explores parameter space by comparing the likelihood of a 
model defined by one set of parameters with that defined by a perturbed set 
of parameters. The likelihood L of a model is defined as 


L = e c 1 jC = --J2^ 1 ^ 
1 2 ^ 


crt 


( 9 ) 


where Xi is the ith data point, /i, is the ith model point, and of is the variance 
of the ith data point, calculated from Table 1 of IE11. Models with high L 
values correspond to better fits to the data than models with low L values, 
and thus the sets of parameters that define the better-fit models are closer to 
the true values of these parameters than the sets of parameters that define 
the worse-fit models. 

At every link of the MCMC, one of the four parameters is chosen randomly 
to be perturbed by adding onto the original parameter value a perturbation, 
the magnitude of which is drawn from normal distributions of certain stan¬ 
dard deviations: 10 4 kg for M 0 ; 0.5 /mi for r m and r 0 , and 1 for /. The 
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likelihood of this perturbed set of parameters is then calculated according 
to Eq. [9j If the likelihood of the model defined by the perturbed set is 
greater than that of the original set, then the perturbed set of parameters 
are accepted as a new link in the chain and the new standard with which 
subsequent sets of perturbed parameters are compared; otherwise, the per¬ 
turbed set is accepted if the ratio of the likelihood of the model defined by 
the perturbed set to that of the original set is greater than a random number 
drawn from a uniform distribution from 0 to 1. In other words, if the L value 
of the model defined by the perturbed set is very low, i.e. a very bad fit to 
the data and thus a very low probability that the data is described by this 
model, then it is highly unlikely that it will be accepted. By calculating the 
likelihood of models defined by different sets of parameters with a bias pro¬ 
portional to L, the MCMC method samples the probability density function 
(PDF) of the parameters. The PDF of each parameter is then the histogram 
of the accepted parameter values. 

We use a total of 2 x 10 6 links in the MCMC, and the parameter PDFs 
are constructed from taking the histograms of the accepted values. Sensi¬ 
tivity tests are conducted by halving and doubling the widths of the normal 
distributions from which the perturbations at each step of the MCMC were 
drawn, with no major changes to our results. See Appendix A for further 
details on these sensitivity tests. 

We apply loose bounds to the parameters: M 0 is restricted to positive 
values; r 0 is bound between 0.2 and 20 /mi and r m is limited to values between 
0.08 and 14 /irn to ensure that the radius grid captures most of the size 
distribution; and / is constrained between 0 and 10 to avoid size distributions 
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narrower than the resolution of the radius grid. The best fit parameters of 
IE11 are well within these bounds. If there arises a perturbation during a 
step of the MCMC that resulted in values of the parameters outside their 
bounds, then the perturbation is redrawn from its normal distribution until 
a value within the bounds results. 

5. Results 

Figure [2] shows the results of the MCMC calculations in the form of scatter 
plots of the accepted sets of parameters for each parameter pair. Areas of 
parameter space with more accepted values (points) correspond to places of 
higher likelihood, and thus better fits to data. It is immediately clear that 
there are multiple families of solutions, which have been assigned different 
colored points. For example, the M 0 vs. r m plot shows at least seven local 
maxima of accepted parameter sets, though the local maximum with the 
lowest M 0 appears to be separate from the others. Conversely, the M 0 vs. 
ro and Mo vs. / plots show only three local maxima at distinct locations 
in parameter space, one at low Mo and two at high Mo. These two high 
M 0 families of solutions can be further split in terms of r 0 and /, with one 
at high r 0 and high / (and slightly lower M 0 ) and one at low r 0 and low / 
(and slightly higher Mo). From the r m vs. ro and r m vs. / plots we see 
that the local maximum with the second lowest r m is resticted to high ro 
and high / values, and that it is the only local maximum out of the seven 
to exhibit this behavior. Therefore, the seven local maxima of r m can be 
distributed among the three local maxima in the other plots as follows: the 
local maximum with r m < 0.8 /rm defines the low M 0 family of solutions 
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(red); the local maximum with 0.8 /im < r m < 1.25 [irri defines the high 
Mo, ro, and / family of solutions (blue); and the remaining local maxima 
with r m >1.25 /irri defines the high M 0 and low r 0 and / family of solutions 
(green). In other words, the red points are solutions corresponding to low 
particulate mass plumes composed primarily of small aggregates with N m > 
100; the blue points are high particulate mass plume solutions consisting of 
large aggregates made up of large monomers, such that they also have N m > 
100; and the green points are also high particulate mass plume solutions, but 
they are comprised of a mixture of spheres and large aggregates composed 
of very few and very large monomers, which we term the “sphere-aggregate” 
solutions. 

Despite the distinctiveness of these solutions in parameter space, most 
of their x m values are >1.5, which is outside the validated range of the T08 
model. Therefore, caution must be exercised when interpreting these re¬ 
sults. As there are currently no quantitaive estimates of the divergence of 
this approximate model from the detailed T-matrix calculations outside the 
validated range (Mishchenko et ah, 1996), we compare our results to that 


of the fractal aggregate scattering model of Rannou et ah (1999, henceforth 
R99). The R99 model is semi-empirical and derives particle scattering prop¬ 


erties using interpolated values from the mean field approximation of Botet 


et al. (1997), which itself is based on an exact theory of scattering by ag¬ 


gregate particles (Xu, 1995). The R99 model has been validated to x m ~ 8 
(Coustenis et al. 2001 Rannou et al. 2002, 2004, 2010; Lavvas et al., 2010[ 
for e.g.) and shows good agreement with the T08 model within the latter’s 
validated range. Thus, it can be used as an approximate fiducial to evaluate 
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the accuracy of the T08 model beyond this range. Only rough agreement is 
necessary however, as the R99 model is more error-prone than the T08 model 
due to nonlinear effects in light scattering (T08). 

Figure [3] shows the percent difference in the scattering efficiency Q sca (top) 
and the phase function P(6>) (bottom) (see Eq. [6| between the R99 and T08 
models for relevant scattering angles and size parameters. The comparison is 
done in the CLR wavelength channel with 300 monomers assumed for each 
case, which is expected for the small and large aggregate solutions given their 
r m and r 0 values; the refractive indices used are those of water ice (Section 
|3| and the fractal dimension is set to 2. The x m values of each of the colored 
phase curves in the bottom panel are indicated by the same colored points 
in the top panel. The comparison shows that the Q sca and P(0) values for 
the two models agree to within 20% for most x m < 3.5, thereby extending 
the range of validity of the T08 model. For the small aggregates, x m ~ 3 for 
the CLR channel (assuming r m ~ 0.3 /rm), and thus it is a valid family of 
solutions. Conversely, the large aggregates have x m ~ 10 for the CLR channel 
(assuming r m ~ 1 /rm), which is well outside the extended validity range, 
and therefore calls into question the existence of this family of solutions. An 
additional issue with the large aggregate solutions is its high r 0 value (~ 
20 /rm), which goes against previous estimates of the sizes of the Enceladus 


plume particles of a few /rm (Kieffer et al. 2009 Hedman et al. 2009). The 
sphere-aggregate solutions are largely independent of the limitations of the 
T08 model, as most of the forward scattering is due to spherical particles. 

From this exercise it is clear that the small and sphere-aggregate solutions 
are likely real, but that the large aggregate solutions may not be. Therefore, 
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we will ignore the large aggregates in the remainder of this work. 

The left side of Figure [4] shows the PDFs of the four varied parameters for 
the small aggregate solution, generated by marginalizing the 2D distributions 
of Figure [2] along one dimension for the small aggregate solution only and 
normalizing the area of the resulting ID histograms. The most probable 
values - the maximum likelihood estimates - correspond to the peaks of the 
PDFs, marked by the green lines. Models generated from these parameter 
values are most likely to describe the data. For the M 0 , r m , and r 0 PDFs, the 
68% (la) confidence intervals are enclosed by the blue lines, which indicate 
the smallest interval in parameter space that encloses 68% of the accepted 
values. The PDF of / does not have a well-defined peak, and thus a lower 
limit is given, with the 68% confidence interval now defined as the 68% 
of accepted values immediately below the most probable value. Table [l] 
gives these most probable values and the associated 68% confidence intervals, 
which we set to be the ranges of uncertainty in our fit to data. 

The small aggregate plume solution yields a most probable ro consistent 
with the r 0 of the spherical particles of IE11, but the total plume mass is 
~ 6 times lower. Meanwhile, the PDF of / is very broad, meaning that 
this parameter is not well-constrained by the Cassini ISS data. The lack of 
a upper limit allows for large / values, leading to narrow size distributions 
centered at r 0 , with low abundances of both smaller and larger particles. 
This is consistent with a scenario where the smaller monomers are depleted 
to form small aggregates of median radius ro- 

The right side of Figure [4] shows the PDFs of the sphere-aggregate so¬ 
lutions. As before, the 68% confidence intervals are shown for PDFs with 
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well-defined peaks (M 0 and /), while for the r 0 PDF only the lower limits 
are given. The r m PDF is multimodal and thus does not have a well-defined 
most-probable value or confidence interval. This is likely due to the lim¬ 
itations of the T08 model, as the corresponding x m values are far greater 
than even the extended validity range given above. Beyond r m > 6 /rm, 
spherical particles dominate the particle population, and thus r m becomes 
unconstrained. 

The results of IE11 indicate that, for spherical particles, mass distri¬ 
butions with / = 1 to 2 best fit the data. This is consistent with our 
sphere-aggregate particle plume results. Our M 0 and r 0 estimates for the 
sphere-aggregate solutions are also consistent with IE11, though the PDF 
for ro has a significant tail towards larger values. This highlights a deficiency 
in forward scattering measurements, as more massive particles are difficult 
to observe due to their low surface to volume ratio. Therefore, the impact of 
massive particles on the phase curves is minimal, allowing the results to be 
poorly constrained at large particle sizes. 

Figure [5] gives representative best fits to the VIO (top), CLR (middle), 
and IR3 (bottom) channels of the Cassini ISS data for the small aggregate 
(red) and sphere-aggregate (green) particle solutions, the latter of which 
reproduces the lower left panel of Figure 6 of IE11. In our case, the small 
aggregate solution provides the best fit, as defined by its reduced chi-square 
value 


where N 


Xr = 


E 


(xi - HiY 


N ~ n ~ 1 T 

10 is the number of data points being fitted to, n 


( 10 ) 

4 is the 
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number of parameters, and the variables in the summation are the same as 
that of Eq. [9j Lower Xr values indicate better fits, though if Xr ^ 1 then 
the error bars on the observations are likely overestimated. Xr = 0-71 for the 
small aggregate solution fit and Xr — 1-5 for the sphere-aggregate solution 
fit. Both solutions give similar fits for the VIO and CLR channels, but 
differ significantly in the IR3 channel at scattering angles < 2°. Whereas the 
small aggregate model curve features a broad forward scattering peak, the 
sphere-aggregate model curve’s forward scattering peak is much narrower 
in IR3. The width of the forward scattering peak is related to tq (IE11), 
but also r m in the case of aggregates (Lavvas et ah, 2010); large particles 
result in narrower/sharper peaks while smaller particles result in broader 
peaks. Therefore, as the small aggregate particles have small r m values, 
their forward scattering peaks are the most broad, increasing R{9) at larger 
scattering angles and prompting low M 0 values. Meanwhile, the lack of small 
monomers in the sphere-aggregate particles narrows the forward scattering 
peak, resulting in low values at large 9, and a higher Mo than for the small 
aggregate solution. 

It is interesting to note that only the IR3 channel offers a way to discrim¬ 
inate between the two solutions. This is caused by the strong dependence 
of the intensity of scattered light on the particle size parameter x m . In our 
fits there are two ranges of particle sizes to consider: the small aggregate 
monomer radius ~ 0.3 p m and the median particle radius for the small ag¬ 
gregate and sphere-aggregate solutions ~ 4 pm; the intensity of the scattered 
light will depend on a combination of the particle radius and the monomer 
radius, if applicable. For the VIO and CLR channels (A = 420 and 635 
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nm, respectively), both the small aggregate and sphere-aggregate particles 
are much larger than the wavelengths considered, leading to narrow forward 
scattering peaks. However, for the IR3 channel (A = 918 nm), the small 
aggregate monomer radius is now much smaller than the wavelength, lead¬ 
ing to a much broader forward scattering peak. Similarly, the low / value 
of the sphere-aggregate solution means that there exists a great number of 
small particles that are responsible for much of the scattering, which are now 
much smaller than the channel wavelength, again leading to a wider forward 
scattering peak. 

In obtaining two separate families of solutions, we have shown that de¬ 
generacy exists in the ISS data from IE11 when both spherical and aggregate 
particles are considered. This degeneracy can be broken by an independent 
mass measurement, such as from Cassini CDA. 

6. Discussion 

We can calculate the solid to vapor mass ratio of the plume for the small 
aggregate solution following the procedure described in IE11. We derive a 
value for the particulate column mass abundance M p ol by relating it to the 
I/F values of the Cassini ISS NAC image from IE11 using their Eq. 8, 

T M p 

F = lZ m = Ko{9)M ^ (11) 

where R(9) is defined in Eq. [7j Dividing M p ol by the vapor column mass 
abundance, M^ ol , at the same altitude above Enceladus then gives the solid 
to vapor mass ratio of the plume at that altitude. From the small aggregate 
curve (red) in Figure [HJ we interpolate a R(6) value of 6.4 x 10 3 km 2 for the 
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CLR channel at 9 = 2.38°, the wavelength channel and scattering angle of the 
NAC image, respectively. This is slightly higher than the value given in IE11, 
5.4 x 10 3 km 2 , owing to the different model curves used for the interpolation. 
Setting M 0 to (25 ± 4) x 10 3 kg for the small aggregate solution and I/F to 
0.07, the brightest pixel value of the NAC image above background (as with 
IE11), we find that M^ ol = (2.7 ± 0.4) x 1CT 7 kg m -2 for that pixel, which 
is ~ 7 times lower than the M/! ol value calculated by IE11. The altitude 
above Enceladus that corresponds to that pixel is similar to that probed by 


the occultations carried out by Hansen et al. (2011) to measure M/ ol using 
Cassini UVIS. Therefore, the solid to vapor mass ratio can be obtained by 
dividing the M// ol value derived here by the M/ ol value derived from averaging 


the second column of Table 1 of Hansen et al. (2011). This gives a solid to 
vapor mass ratio for the small aggregate plumes of 


Ko, 

Ml, 


0.07 ±0.01, 


( 12 ) 


which is again 7 times lower than that of IE11. The error arises from uncer¬ 
tainties in both the solid and vapor column mass abundances and is likely 
underestimated, as we have assumed zero uncertainty for the R(9 = 2.38°) 
and I/F values. It should be noted that, while K 0 (9 ) is defined using all 
available WAC images, the comparison with UVIS requires that we use just 
the brightest pixel of the one NAC image to set the I/F value. In other 
words, the global I/F values of the plume in the WAC images are assumed 
to be the same as that of the brightest plume pixel in the NAC image, which 
is equivalent to assuming that the particle size distribution is independent of 
altitude above Enceladus. This is only an approximation, as evidence exists 
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that the particle size distribution does change with altitude (Hedman et al. 


2009 Postberg et ah, 2011). However, for simplicity we do not consider this 


effect in this work. 


Our calculated solid to vapor ratio is consistent with that of Kieffer et ah 


(2009), who provided an upper bound of ~ 0.1 — 0.2 using observations from 


Cassini ISS and UVIS (Porco et al. 2006). It is also consistent with the 
lower end of the range of solid to vapor ratios (~ 0.01 — 100) calculated by 


Hedman et ah (2009) using Cassini VIMS observations and spherical and 


irregular particle scattering models, though aggregates were not considered. 


The K 0 (6 ) quantity in Eq. 11 is extremely useful as it can readily convert 
an observed I/F value into the column mass abundance. It also varies de¬ 
pending on the particle scattering properties, wavelength of scattered light, 
and particle size distribution. Figure [6] shows the K o (0) values of the small 
(solid lines) and sphere-aggregate (dashed lines) solutions at scattering an¬ 
gles between 0° and 50° for the VIO (blue), CLR (green), and IR3 (red) 
wavelength channels. The small aggregate K o (0) values are, as expected, 
higher than that of the sphere-aggregate solution, since there is less mass for 
the same intensity of forward scattering. The sphere-aggregate K o (0) val¬ 
ues are typically within 10% of those of IE11, further reinforcing the lack of 
dependence of the sphere-aggregate solutions on the limitations of the T08 
model. A full set of K o (0) values for the two families of solutions and the 
best fit solution of IE11 at the Cassini WAC and NAC channel wavelengths 
can be found in the online supplemental material. 

Aside from what can be retrieved from observations, it is also impera¬ 
tive that the particle solutions are physical and adhere to theoretical con- 
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straints of formation and evolution. As IE11 has already considered high 
mass plumes, we will focus on the low mass small aggregate solutions and 
their possible formation mechanisms for the remainder of this work. 

Aggregate particles in the Solar System typically form through the co¬ 
agulation of monomers, a process that depends heavily on the free-floating 


monomer number density (Lavvas et al. 2010, for e.g.,), i.e. those not already 
incorporated into aggregates. Therefore, much of the coagulation process 
must take place within the plume vents, where the particle number density 
is higher. This is consistent with our small aggregate particle solution, which 
shows a low abundance of monomers compared to aggregates in the plume, as 
the monomers must have already coagulated to form aggregates before they 
were ejected into space. The residence time of particles inside the plume 
vents is of order D/v, where D is the depth of the liquid vapor interface and 
v is the average velocity of the particles. In order for complete coagulation 
of monomers within the plume vents, this residence time must be greater or 
equal to the coagulation time scale, approximated by 


1coag 


(1 

dn 

\ n 

dt 


-1 


1 

nK’ 


(13) 


where n is the free-floating monomer number density in the vent, and K is 
the coagulation kernel. In the free molecular limit (Knudsen number >> 1), 
the Brownian coagulation kernel for two spherical particles with radii rq and 
r 2 is 


K = (r 1 +r 2 y 


6kT 


\ Pp 


1 1 
+ 


r 3 

Jl 


r 3 
' 2 , 


(14) 


where k is the Boltzmann constant, T is temperature, and p p is the solid 
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particle density (0.917 g cm 3 ). Setting rq = r 2 = r m , 


K = 


\ 


192 kTr n 
Pp 


( 15 ) 


(Pruppacher and Klett, 1978). To find a value for the coagulation time scale, 


we consider the rate at which particulate mass is emerging from the vents, 
which can be expressed as 


= \~ T t,Pt nAv = PmAv, 


(16) 


where p m is the total mass of particulates per unit volume and A is the total 
vent area of the Enceladus plumes. 1E11 gives this rate as dM p /dt = 51 kg 
s -1 , but the reduced plume mass of the small aggregate solutions results in a 
smaller value. Rather than assuming that dM p /dt is linearly proportional to 
M 0 however, we instead appeal to the solid to vapor mass ratio by considering 
the rate at which vapor is emerging from the vents, dM v /dt , which has been 


determined by Hansen et al. (2011) to be ~ 200 kg s 1 based on Cassini 


UVIS observations. Thus, 


dM v 

dt 


= p v Av = 200 kg s 


-i 


(17) 


where p v is the density of the vapor at the liquid-vapor interface, 4.85 x 
10~ 3 kg m -3 at T = 273 K. Dividing Eq. 16 by Eq. 17 gives us the solid 


to vapor mass ratio ~ [dM p /dt / dM v /dt] and allows us to eliminate A and 


v. Substituting the result into Eq. 13 then allows us to eliminate n. The 
necessary condition for the completion of monomer coagulation before exiting 
the vent then becomes 
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D 

v 


1 /4 


> T coag — ( 0 7Tr m 


p p dM v jdt 
p v dM p /dt 


20 s 


(18) 


where we use r m = 0.3 pm and T = 273 K. For [dM p /dt / dM v /dt] we use 
the solid to vapor mass ratio derived in Eq. 


Several assumptions have gone into Eq. [18j One is that the particles and 
the vapor emanate from the same liquid-vapor interface, or at least very close 
above it in the vent. This is necessary to maintain high T and p p values. If 
the particles form in flight higher up in the vent, as in a convective cloud, 
the reduction in p p due to lower T and/or condensation on the vent walls 
would drastically increase T coag . The fact that some of the plume particles 


are Na-rich (Postberg et al. 2009, 2011) argues that they originate from 
the liquid, but we cannot rule out a vapor origin for the Na-poor particles. 
Another assumption is that the different parts of the erupting area A behave 
the same way whether they are a single crack running the length of each 


tiger stripe, a series of 100 discrete sources (Porco et ah, 2014), or some 


combination of the two (Tian et al. 2007 Hansen et al. 2011; Postberg 


et al., 2011; Spitale et al., 2015, for e.g.). A third assumption is that the 


solid to vapor mass ratio remains the same as the plumes vary with orbitial 


phase (Hedman et al. 2013). These assumptions are difficult to test, but 
they seem reasonable. 

A fourth assumption is that the relevant vapor density and velocity are 
close to those at the liquid-vapor interface and not those at the top of the 
vent where the particles and vapor exit to space. From the vertical distri¬ 
bution of particles exiting the vent, one infers velocities of order 60-90 m 
s -1 (Porco et ah, 2006; IE11). However, the particles are accelerated to this 
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speed within a few meters of the surface (Schmidt et al., 2008 Ingersoll and 


Pankine 2010), so the exit speed underestimates the time spent in the vent 
during which coagulation takes place. Postberg et al. (2009) use 500 m s' 1 
for the speed of the vapor and the triple point value for its density, but 


Nakajima and Ingersoll (Submitted) point out that these are overestimates. 


Their model includes friction with the walls of the crack as the gas flows 
upward. This produces backpressure at the liquid-vapor interface, and for a 
long, narrow crack the backpressure severely limits the evaporation rate and 
the upward velocity. For crack widths less than 0.1 m, the upward velocity 
is less than 10 m s _1 below depths of 1 km. If the evaporating surface is 
another 1 km below that level, the time spent in the vent is well over 100 s 
and the condition for coagulation of aggregates (Eq. [18]) is comfortably sat¬ 
isfied. Crack width is an important parameter, since it affects the velocity. 


Nakajima and Ingersoll (Submitted) derive their estimate by matching the 


ratio of the power released through latent heat to that released through in¬ 
frared radiation, but one should remember that these numbers are uncertain 
by factors of 2 or more. 

Our results favor a less violent source for the plumes than the earlier 
results of IE11. Here, “less violent” means steady, controlled evaporation 
rather than explosive boiling. For example, it could be that most of the ice 
particles condense directly from the expanding vapor above the liquid-vapor 
interface rather than erupt as a bubbly liquid, which then breaks up and 
freezes into small ice particles in the vacuum of space. Alternatively, the 
particles could be generated by bubbles of methane or CO 2 rising through 
liquid water and breaking at the surface, thereby sending up a fine spray of ice 
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particles. Porco et al. (2006) mentioned two extremes: particles condensing 
directly from the vapor and particles forming from the breakup of a boiling 
liquid, the latter analogous to a cold Yellowstone geyser. The relatively large 
solid to vapor ratio of IE11 favors the violent source, while the lower ratio of 
our small aggregate solution allows for less violent sources. Thus, the present 
results make the ISS data compatible with a wider variety of sources than 
the IE11 results. 


Schmidt et al. (2008) and Ingersoll and Pankine (2010) used hydrody¬ 


namic models to simulate the condensation of ice particles from water vapor, 
and obtained solid to vapor mass ratios of 0.05 - 0.06 and 0.015, respec¬ 


tively. The lower value calculated by Ingersoll and Pankine (2010) results 


from the condensation of vapor onto the vent walls, which Schmidt et al. 


(2008) did not consider. These values of the solid to vapor mass ratios are 
more consistent with our small aggregate solution than with the results of 


IE11. Kieffer et al. (2006) considered explosive decomposition of a hydrate 
clathrate that contains methane and other gases in much greater abundances 
than their solubility in liquid water would allow. The solid to vapor mass 
ratio of such an event is uncertain, but it could be lower than that quoted by 
IE11 and closer to the ratio derived here from the small aggregate solution. 


Meanwhile, Hsu et al. (2015) and Postberg et al. (2009, 2011) showed that 
silica grains and sodium salts are present in both E ring and plume parti¬ 
cles, which argues strongly for frozen droplets originating from a salty liquid 


rather than ice particles forming directly from vapor condensation. Postberg 


et al. (2009) suggested that ascending bubbles of plume gases (C0 2 , N 2 , CO, 


CH 4 ) can disperse the liquid droplets into the vapor, and that the evapora- 
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tion of the boiling liquid source is taking place over a large horizontal area 


so that freezing of the vents is suppressed. Nakajima and Ingersoll (Submit¬ 


ted) showed that the backpressure due to friction on the walls of the channel 


could severely limit the evaporation rate. Ingersoll and Nakajima (Submit¬ 


ted) argue that this “controlled boiling” allows bubbles of vapor to break at 


the surface instead of throwing up large amounts of spray as in boiling into 
vacuum. All of this points to a relatively gentle source and a relatively low 
solid to vapor mass ratio. 

Our small aggregate solution does not rule out a vapor-based source for 


the plume particles. Parkinson et al. (2008) showed that spherical particles 
with radii similar to r m can form from water vapor nucleating on involatile 
grains in a few seconds, while spherical particles with radii up to 3 /mi can 
form by condensation in about 100 s, though large condensation nuclei and/or 
a high density of water vapor may be necessary. Similarly, non-spherical ice 
crystals resembling aggregates, e.g. dendritic shapes, can form directly from 
condensation and nucleation of water vapor, with growth rates comparable 


to that of spherical particles (Pruppacher and Klett, 1978). Bullet rosettes 
- highly non-spherical, polycrystalline ice particles made up of clumps of 
hexagonal columns - can also yield high scattering cross sections for a small 
particle mass, but they are formed by the fracturing of ice during the freezing 


of a large water drop (Baum et al. 2011). 

We have shown that forward scattering observations of the Enceladus 
plumes from Cassini ISS can be well-fit by plumes made up of ice aggregates 
and solid ice spheres. This results in a bifurcation in the allowed total partic¬ 
ulate mass of the plumes: small aggregate plumes are six times less massive 
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than plumes made of larger aggregates and spheres. The small aggregate 
plumes, in particular, lead to a solid to vapor mass ratio of 0.07 ± 0.01, 
which is suggestive of a “gentle source” for the plume particles. 
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Figure 1: A schematic of the plume particle mass distribution dM/dr as a function of the 
particle radius r, with the median particle radius ro and minimum aggregate radius r m i n 
labelled. 
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Figure 2: Scatter plots of the accepted sets of parameters from the MCMC calculations for 
each parameter pair. The regions with higher concentrations of points correspond to areas 
of higher likelihood. Red points correspond to the small aggregate family of solutions 
defined by r m < 0.8 p,?n, with low M 0 ; blue points correspond to the large aggregate 
family of solutions defined by 0.8 ^m < r m < 1.25 jim , with high Mq, ro, and /; and 
green points correspond to the sphere-aggregate family of solutions defined by r m > 1.25 
/urn, with high Mq and low Tq and /. 
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Figure 3: Percent difference in the scattering efficiency Q sca (top) and the phase function 
P(0) (bottom) between the models of Rannou et al. ( 1999| ) and |Tomasko et ah (20081 for 
relevant values of the size parameter x m and scattering angle 6. The comparison is done 
in the CLR wavelength channel with 300 monomers assumed for each case. The refractive 
indices used are those of water ice (Section [3]) and the fractal dimension is set to 2. The 
x m values of each of the colored phase curves in the bottom panel are indicated by the 
same colored points in the top panel. The gray shaded region in the top panel indicates 


the range in x m for which the model of Tomasko et al. (2008) has been validated. 
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Figure 4: The probability density functions (PDFs) of the total mass of the plume 
Mq (top), the monomer radius r m (middle-top), the median particle radius r o (middle- 
bottom), and the width of the size distribution / (bottom) for the small aggregate (left) 
and sphere-aggregate (right) particle plume solutions. Note the different abscissa val¬ 
ues between the rg and r m plots. Each histogram contains 50 bins. The most probable 
value is marked by the green line. For the Mo, small aggregate r m and ro, and sphere- 
aggregate / PDF’s, the 68% confidence intervals are enclosed by the blue lines, while for 
the small aggregate / and sphere-aggregate r 0 PDFs only the lower 68% confidence in¬ 
terval is marked. The sphere-aggregate r m PDF is multimodal and thus does not have a 
well-defined most-probable value or confidence interval. 


37 













































Figure 5: Representative best fits to the VIO (top), CLR (middle), and IR3 (bottom) wave¬ 
length channel data from Cassini ISS for the small aggregate (red) and sphere-aggregate 
(green) plume particle solutions. Parameters used for the small aggregate solution fit are: 
M 0 = 22.58 x 10 3 kg, r m = 0.331 /im, ro = 3.9 /mi, and / = 7.79. Parameters used for 
the sphere-aggregate solution fit are: M 0 = 172.42 x 10 3 kg, r m = 4.87 /tm, r 0 = 4.81 
/rm, and / = 1.71. 
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Figure 6: Variations of Ko((9) = R(0)/Mo as a function of scattering angle for the VIO 
(blue), CLR (green), and IR3 (red) wavelength channels and the small aggregate (solid 
lines) and sphere-aggregate (dashed lines) solutions. 


39 










Parameter 

Small Aggregates 

Sphere-Aggregates 

M 0 (10 3 kg) 

25 (22 - 30) 

166 (119 - 205) 

r m (pm) 

0.27 (0.24 - 0.35) 

- 

f o (pm) 

4.20 (3.11 - 10.51) 

3.90 (3.34 - oo) 

/ 

2.5 (2.1 - oo) 

1.34 (1.05 - 1.89) 


Table 1: Most probable values and (in brackets) 68% confidence intervals of the retrieved 
parameters, where Mq = total mass of the plume, r m = monomer radius, ro = median 
particle radius, and / = width of the size distribution for the small aggregate (left) and 
sphere-aggregate (right) particle plume solutions. The sphere-aggregate ro and the small 
aggregate / PDFs have unconstrained upper bounds and so only the lower bounds are 
given, defined as the lower limit of the 68% of the accepted values immediately below the 
most probable value. The monomer radius for the sphere-aggregate solution is omitted, as 
it is multimodal and thus does not have a well-defined nrost-probable value or confidence 
interval. 
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Appendix A. Sensitivity Tests 


A major uncertainty in the Markov Chain Monte Carlo method is the 
standard deviations of the normal distributions from which the perturba¬ 
tions at each step are drawn. For our nominal case we use 10 4 kg for Mo, 0.5 
fin 1 for r m and tq, and 1 for /, as they are on or smaller than the value of 
the actual quantity itself. To test the sensitivity of our results to the stan¬ 
dard deviations, MCMC runs are conducted where these values are halved 
or doubled for all parameters, while leaving all other properties of the runs 
unchanged. 

Figures [AT] and [M2] show our results for the small aggregates and sphere- 
aggregates, respectively. In comparison with the nominal case (gray shaded 
regions), the halved and doubled cases are different but maintain much of 
the same shapes and widths. Table A.l presents the most probable values 
and 68% confidence intervals of the optimized parameters for the halved and 
doubled cases, which are in satisfactory agreement with that of the nominal 
case (Table [I]) to within its own 68% confidence intervals. Therefore, we can 
conclude that our nominal results are robust. 
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Half 


Double 



Figure A.l: Same as the left side of Figure [4] (the small aggregates), but for the sensi¬ 
tivity tests where we half (left) and double (right) the standard deviations of the normal 
distributions from which we draw the perturbations to the parameters at each step of the 
MCMC. The most probable values and 68% confidence intervals are given by the green 
and blue lines, respectively. The parameter PDFs shown on the left side of Figure [4] are 
plotted in gray for comparison. 
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Double 






/ / 

Figure A.2: Same as Figure |A.1[ but for the sphere-aggregates. The parameter PDFs 
shown on right side of Figure [4] are plotted in gray for comparison. 
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Parameter 

Small Aggregates 

Sphere-Aggregates 

Half: 

M 0 (10 3 kg) 

29 (23 - 31) 

169 (125 - 230) 

r m (pm) 

0.26 (0.23 - 0.34) 

- 

f o (pm) 

4.07 (3.19 - 13.00) 

4.24 (3.51 - oo) 

/ 

2.9 (2.2 - oo) 

1.37 (1.05 - 1.77) 

Double: 



M 0 (10 3 kg) 

28 (23 - 31) 

170 (118 - 247) 

r m (pm) 

0.26 (0.23 - 0.34) 

- 

r 0 (pm) 

4.17 (3.25 - 12.70) 

4.67 (3.67 - oo) 

/ 

2.8 (2.2 - oo) 

1.36 (1.08 - 1.78) 


Table A.l: Same as Table [I] but for the sensitivity tests where we half (top) and double 
(bottom) the standard deviations of the normal distributions from which we draw the 
perturbations to the parameters at each step of the MCMC. 
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